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In this theoretical paper, we investigate coherence properties of the near-resonant light scattered 
by two atoms exposed to a strong monochromatic field. To properly incorporate saturation effects, 
we use a quantum Langevin approach. In contrast to the standard optical Bloch equations, this 
method naturally provides the inelastic spectrum of the radiated light induced by the quantum 
electromagnetic vacuum fluctuations. However, to get the right spectral properties of the scattered 
light, it is essential to correctly describe the statistical properties of these vacuum fluctuations. Be- 
cause of the presence of the two atoms, these statistical properties are not Gaussian : (i) the spatial 
two-points correlation function displays a speckle-like behavior and (ii) the three-points correlation 
function does not vanish. We also explain how to incorporate in a simple way propagation with a 
frequency-dependent scattering mean-free path, meaning that the two atoms are embedded in an av- 
erage scattering dispersive medium. Finally we show that saturation-induced nonlinearities strongly 
■ modify the atomic scattering properties and, as a consequence, provide a source of decoherence in 

multiple scattering. This is exemplified by considering the coherent backscattering configuration 
where interference effects are blurred by this decoherence mechanism. This leads to a decrease of 
04 ■ the so-called coherent backscattering enhancement factor. 
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I. INTRODUCTION 

> 
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Over the past ten years, cold atomic gases have gradually become a widely employed and highly tunable tool for 
testing new ideas in many areas of quantum physics: quantum phase transitions (Bose-Einstein condensation, Fermi 
degenerate gases, Mott- Hubbard transition) jHQ.13- quantum chaos applications in metrology Q, disordered 
systems 0, Q to cite a few. In the latter case, cold atomic vapors act as dilute gases of randomly distributed atoms 
multiply scattering an incident monochromatic laser light. In this case, the scattered light field exhibit a speckle-like 
structure due to (multiple) interference between all possible scattering paths. The key point is that the disorder 
1-^ I average is insufficient to erase all interference effects. This gives rise to weak or strong localization effects in light 
transport depending on the strength of disorder 0, . A hallmark of this coherent transport regime is the coherent 
backscattering (CBS) phenomenon: the average intensity multiply scattered off an optically thick sample is up to 
twice larger than the average background in a small angular range around the direction of backscattering, opposite 
to the incoming light H3. This interference enhancement of the diffuse reflection off the sample is a manifestation 
q-i of a two- wave interference. As such, it probes the coherence properties of the outgoing light [llj . The CBS effect in 
cold atomic gases has been the subject of extensive studies in the weak localization regime, both from theoretical and 
. ! experimental points of view |l2j . In particular, modifications brought by atoms, as compared to classical scatterers, 
for light transport properties (mean-free path, coherence length, CBS enhancement factor) have been highlighted. 
?H I They are essentially due to the quantum internal atomic structure 0, Q] . 

Another interesting feature of atoms is their ability to display a nonlinear behavior: the scattered light is no more 
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wide variety 01 p 

self-focusing effects, dynamical instabilities, etc flol llH IT7I Il8|. For a weak nonlinearity, introducing an intensity- 



proportional to the incident one. This leads to a wide variety of phenomena, like pattern formation, four- wave mixing, 

ing an inten 

dependent susceptibility is enough to properly describe these effects, including quantum properties [TBI llH , l2(i | . e.g. 
the Kerr effect (intensity dependence of the refractive index) can be obtained with a \^ nonlinearity. However, when 
the incident intensity is large enough, and this is easily achieved with atoms, perturbation theories eventually fail 
and a full nonlinear treatment is required. For a single two-level atom, the solution is usually given by the so-called 
optical Bloch (OB) equations. Together with the qua ntum regression theorem, they allow for a complete description 
of the spectral properties of the fluorescence light 22]. In particular, these equations show that the atomic nonlinear 
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behavior is intrinsically linked to the quantum nature of the electromagnetic field. More specifically, as opposed 
to classical nonlinear scatterers, the radiated light exhibits quantum fluctuations characterized by peculiar time 
correlation properties. They define a power s pect rum, known as the Mollow triplet, emphasizing inelastic scattering 
processes at work in the emission process |22l l23l 124) . 

However, even if all these aspects are well understood in the case of a single atom exposed to a strong monochromatic 
field |22j , the situation changes dramatically in the case of a large number of atoms where a detailed analysis including 
both quantum nonlinear properties and coherence effects is still lacking. Until now, the nonlinear coupling between the 
atoms and the quantum vacuum fluctuations is either included in a perturbative scheme or simply described 

by a classical noise [2^, |2(| |27l IH, l^j . In the dilute regime A <C R where the light wavelength A is much less than 
the average particle separation R, one expects the quantum fluctuations to reduce the degree of coherence of the 
scattered light. This will alter not only propagation parameters (mean- free path, refraction index), but also weak 
localization corrections to transport, and the CBS enhancement factor, which is related to the coherence properties 
of the scattered light field Q, El- We want here to stress that, even beyond interference and weak localization 
phenomena, any transport property which may be influenced by saturating the atomic transition deserves a special 
and necessary study on its own. The most striking systems falling in this category where both nonlinear and disordered 
descriptions are intimately interwoven are coherent random lasers [30| where interference effects lead to localized light 
modes inside the disordered medium, comparable to resonator eigenmodes in standard lasers. Even if, in this case, 
one would require an active {i.e. amplifying) medium, a key point is the understanding of the mutual effects between 
multiple interferences and nonlinear scattering. 

In the present paper, we will focus on the rather simple case of two atoms in vacuum. Our aim is threefold: (i) firstly 
to properly calculate quantum correlations between pairs of atoms as a crucial step towards a better understanding 
of the physical mechanisms at work; (ii) secondly to implement a method allowing for a simple incorporation of 
frequency-dependent propagation effects; (iii) finally to understand, in the CBS situation, the modifications brought 
by the (quantum) nonlinearity to the interference properties. We hope that these points, once mastered, can provide 
an efficient way to produce realistic computer models to simulate real experiments. Point (i) alone could easily be 
solved using the standard OB method pul But the latter almost becomes useless regarding point (ii), since 

frequency-dependent propagation leads to complicated time-correlation functions. From a numerical point of view, it 
also leads to such large linear systems of coupled equations that its practical use is limited up to only a few atoms, 
very far from a real experimental situation. For these reasons, we will rather use the quantum Langevin method 
for our purposes. This method not only solves points (i) and (ii), but also leads to a simple explanation of point 
(iii), through a direct evaluation of the quantum noise spectrum. Note however that, in the absence of any effective 
medium surrounding the two atoms, and as long as only the numerical results are concerned (but not the physical 
interpretation), the quantum Langevin approach is completely equivalent to solving the multi-atoms optical Bloch 
equations like in 0, [^3 . 

This paper divides as follows: in section [B] the notations are defined and the quantum Langevin approach is 
explained for the single atom case. In section ITTT1 the method is adapted to the case where two atoms are weakly 
coupled by the dipole interaction. The validity and relevance of the method is controlled by a comparison with a 
direct calculation using OB equations. Then, in the CBS configuration, numerical results for different values of the 
laser intensity and detuning are presented and discussed. In particular, possible reasons for the reduction of the 
enhancement factor are put forward. 

II. SINGLE TWO-LEVEL ATOM CASE 

A. Time-domain approach 

We consider an atom with a zero angular momentum electronic ground state ( J g = 0) exposed to a monochromatic 
light field. The light field frequency lul is near-resonant with an optical dipole transition connecting this ground state 
to an excited state with angular momentum J e = 1. The angular frequency separation between these two states is 
loq and the natural linewidth of the excited state is T. We will denote hereafter by 6 = lul — luq the laser detuning. 
The ground state is denoted by |0 0) while the excited states are denoted by |I m e ), with m e = —I, 0, 1 the Zeeman 
magnetic quantum number. As we assume no magnetic field to be present throughout this paper, the excited state is 
triply degenerate. 

In the Heisenberg picture, this two-level atom is entirely characterized by the following set of 16 time-dependent 
operators: 

IF = |00}(00| ; IL e msmi = \lm e )(lm' e \ ; £>+ e = |lm e )(00| ; £>~ e = |00)(1 m e | (1) 



The atomic operators obey the completeness constraint 



1 = IF + n e (2) 

where II 9 and IF = ^ m 11^ m are the ground and excited state atomic population operators. 

The full atom- field Hamiltonian Ti. is the sum of the free atom Hamiltonian Ha — hcuoTl e , of the free quantized 
field Hamiltonian Tip — X)k ej_k ^ k a iLe ak£ an( ^ °^ the dipolar interaction V = d • (E^ + Ey) between the atomic 
dipole d, the classical laser field E^ and the quantum electromagnetic vacuum field Ey. Performing the usual 
approximations of quantum optics, i.e. neglecting non-resonant terms (rotating wave approximation) and assuming 
Markov-type correlations between the atomic operators and the vacuum field, one obtains the quantum Langevin 
equations controlling the time evolution of any atomic observable O in the rotating frame |22l |25| : 



dO i i F 

— = i6 L [o, n e ] - - v+pL+(R) - J2P, z?-]nJ-(R) - - (on e + wo) + r£ v+ov- + r (R, t), 

9 9 9 

(3) 

where il q + (resp. fi^ ) are the components of the Rabi frequency of the positive (resp. negative) frequency parts of 
the incident laser beam, i.e. Ml = — d~E where d is the dipole strength. Finally J~o(t) is the Fangevin force depicting 
the effects of the quantum fluctuations of the vacuum electromagnetic field and reads as follows: 

fo(t) = - l -Y,(-mo,v+}n°+ q (R,t) ^n°-(R,t)[o,v-], (4) 

9 9 

where f2 0+ (R, t) is the vacuum Rabi field operator 

flP+(R,t) = ~ y £{u)eaUt*y^~ i{ "~" L){t ~ to) (5) 

k,e±k 

with to an initial time far in the past. From the preceding expression, one can calculate the time correlation functions 
of the vacuum field 33] : 

(-l)^+(R,*),Q°,-(R,i')] =4T6 qq ,f(t-t% (6) 

where /(r) in a function centered around r = 0, whose width r c is much smaller than any characteristic atomic 
timescale (i.e. r c <$i ui^ 1 -C T _1 ) and whose time integral is equal to unity. Thus, hereafter, f(r) will be safely 
replaced by a 5-function: /(r) — * S(t). 

The time evolution for the expectation values is obtained by averaging over the initial density matrix er(to), i.e., 
(0(t)) — Tr(0(t)a(to)). Since the atom and the vacuum field are supposed to be decoupled initially, a (to) is simply 
Cai(io) ® |0}(0| (|0) being the vacuum field state). Because of the normal ordering, one immediately gets: 

(Fo(t)) = 0, (7) 
and the time correlation functions of the Fangevin forces: 

(^ (^(0> = -r^[o(*),i^(t)][oxo.^(0]^(*-* / )- (8) 

The physical picture of the quantum Fangevin approach is to represent quantum fluctuations by a fluctuating force 
acting on the system, in analogy with the usual Brownian motion. Not surprisingly, this leads to a diffusive-like 
behavior of expectation values. More precisely, because of the 8- function in Eq. (jHJl, we can set t! = t for the atomic 
operators and we finally obtain in the stationary regime t 3> to- 

(J r o(t)fo'(t')) = ~D o'S(t-t% (9) 

where D is a matrix of diffusion constants depending only on the stationary values of the atomics operators. The 
stationary hypothesis also results from the fact that these correlation functions only depend on the time difference 
t-t'. 

From this, it is possible to prove that the quantum regression theorem applies |22ll34| . allowing for the calculation of 
two-times correlation functions of the atomic operators and of their expectation values. From their Fourier transforms, 
one can obtain the spectrum of the radiated light. But, for the reasons mentioned in the introduction, we will explain 
how these properties can be obtained in a much simpler way by directly translating the Fangevin equations in the 
Fourier domain I'll. 



B. Frequency-domain approach 



First, because of the constraint J3J), only 15 atomic operators are actually independent. More specifically, we will 
use the following set, denoted by the column vector X: 



IT =- [IT - IFl 



x <^ 



n 



|lm e )(lm e | m e =/= m' e 
|lm e )(00| 
|00)(lm e | 

The Langevin equations for X then formally read as follows: 



V+ 



(10) 



— X(t) =AfX(t)+L + F(t), 



(11) 



where M is a time- independent matrix depending on the laser Rabi frequency il L± , L is a constant vector scaling 
with r and F(t) is a vector characterizing the Langevin forces at work on the atom (for simplicity, we have dropped 
the explicit position dependence). The stationary expectation values are then simply given by: 



(X) 



-M _1 L. 



(12) 



Using Kubo's notations, the Fourier transforms of the different quantities are defined as follows: 

/[A] = J dtj{ty 



/,/)= / ^/[A]e- A *, 



(13) 



leading to the Langevin equations in the frequency domain: 

(-iAl - M) X[A] =2tt5[A]L + F[A]. (14) 

Introducing the Green's function G[A] = (— iAt — M)~ , the solution of the preceding equations simply reads: 

X[A] =G[A](2tt(5[A]L + F[A]). (15) 

Using G[0] — — M _1 and 1)12(1 . this solution separates into a non-fluctuating part X^ [A] and a fluctuating (frequency- 
dependent) part Xjr[A]: 



X L [A] =2tt<5[A](X) 
X F [A] = G[A]F[A] 



(16) 



From the linearity of the Fourier transform, we still have (F[A]) = implying (Xj?[A]) = 0. The time correlation 
functions for the Langevin force components, Eq. JSJ, become: 



(F P [A']F,[A]) =2tt5[A' + A}D 



pg- 



(17) 



where the 2ttS[A' + A] function is a direct consequence of the time-translation invariance, i.e. that we calculate the 
correlation functions in the stationary regime. This implies that the correlation function for the components of X F 
in the frequency domain are: 



((XHA']) p (X F [A]) ? ) = 27n5[A + A'] (GD t G) pq 

where the superscript t means matrix transposition. 

The field radiated at frequency A by the atom at a distance r 3> A (far- field regime) reads as follows: 



(18) 



n+[A] = 



o p ikr 

2 qq q [ 1 kr ' 



(19) 



where we use implicit sum over repeated indices and where V v is the projector onto the plane perpendicular to vector 
r: 



V\ q , = -e q V r e q , =e q (t- ^ e q , = 8 qq , (-1)^-^, (20) 

where the bar denotes complex conjugation and where (r *r) is a dyadic tensor. 

The correlation functions (Vt~, [A']f2+[A]) of the light emitted by the atoms is then proportional to (X>i [A']2?~ [A]) 
and read as follow: 

(ft-,[A']0+[A]> oc (2tt) 2 S[A}S[A'](V+)(V-) + 2ttS[A' + A] ^ Gi>p<(A')G ip (A)D p ' p , (21) 

p'p 

where the index i (resp. i') corresponds to T>~ (resp. T)^,). The non-fluctuating part gives rise to a spectral component 
of the emitted light at exactly the incident laser frequency and is thus naturally called the elastic part. The fluctuating 
part gives rise to the inelastic Mollow triplet spectrum J35|, whose properties (position and width of the peaks) are 
given by the poles of G[A], i.e. by the complex eigenvalues of M. Actually, we simply recover the results of the 
quantum regression theorem, which states that the atomic time correlation functions evolve with the same equations 
than the expectation values (X) = M (X) + L |22L I2H . 



III. TWO-ATOM CASE 
A. Optical Bloch equations 

We now consider two isolated atoms, located at fixed positions Ri and R2. Defining R = R2 — Ri = Ru (with 
R = |R| and u the unit vector joining atom 1 to atom 2), we assume the far-field condition R ^> A to hold. We also 
assume that R is sufficiently small for the light propagation time R/c to be much smaller than any typical atomic 
timescales (r -1 , <J _1 , fi^ 1 ). In this regime, all quantities involving the two atoms are to be computed at the same time 
t. The contribution of the atom-atom dipole interaction in the Langevin equation for any atomic operator O reads: 



dO 
~dt 



dip. 



= *? {([°'^ + ] v " v \ 7 + [°'^ 2+ ] v * v v) tr + [°^~] +v * +v * [°^A) e -kR-} 

(22) 

In the OB equations, the two-atom system is entirely described by the set of 256 operators Xij made of all possible 
products XjX?. The stationary expectation values (X^) are then obtained as solutions of a linear system resembling 
equation fT2"|) . This is the approach used in (^3 , where such optical Bloch equations are solved. 

Since the two atoms are far enough from each other, the electromagnetic field radiated by one atom onto the other 
can be treated as a perturbation with respect to the incident laser field. More precisely, the solutions (Xy) can be 
expanded up to second order in powers of g and g: 

(X^ = (Xy)® + g + g (Xy>® + gg (X^S) + g* ( Xij )^ + f (X^S) (23) 

where the complex coupling constant g is: 

3Texp(ikR) 

J 2 k,R V ' 

In fact, it will be shown below that both terms in g 2 and g 2 give a vanishing contribution to the coherent backscat- 
tcring signal. 

As explained in the introduction, this approach has two drawbacks: (i) the solutions obtained in this way are global 
and, thus, do not provide a simple understanding of the properties of the emitted light; (ii) when the two atoms are 
embedded in a medium whose susceptibility strongly depends on the frequency, the field radiated by one atom onto 
the other at a given time t now depends on the atomic operators of the first atom at earlier times (since retardation 
effects become frequency dependent). Time correlation functions in the dipole interaction then explicitly show up. 

B. Langevin approach 

The Langevin equations for the two sets of atomic operators X Q , with a = 1, 2, read formally: 

X" = M"X Q + L + F" + gT q+ ^ a Vf q ,V p q r + gV p q + Vf q ,T q '"X a , (25) 



where /3 denotes the other atom and where T 9± are 15 x 15 matrices defined by [X t ,Vf] = ±2T? ± X j . Taking the 
Fourier transform of these equations, one gets: 

X Q [A] = G a [A] (2tt<5[A]L + F a [A]) + gG a [A]T q+ Vf q , (x Q ® dJ~) [A] - gG Q [A]7>*,T 9 '- (Z>f+ ® X") [A], (26) 

where <g) is the convolution operator: 

(4®5)[A] = ^- /"/ dAidA 2 J[Ai + A 2 - A]A[Ai]B[A 2 ]. (27) 



2tt . 

Introducing, for simplicity, the following notations: 

X" <0) [A] = G a [A] (2tt<5[A]L + F q [A]) 
G<[A] = G a [A]T«' + 7>^ , (28) 

G Q ^[A] =G Q [A]r 9 '-p^, 
equation lj2oTl becomes: 

X Q [A] = X Q<0, [A] + gG a i [A] (X Q ®D^") [A] — gG a ? [A] ® X Q ) [A], (29) 
from which one gets the expansion in power of g and g (up to gg) for the atomic operators: 

. ^!r.,^, n W „^J°)\,., „ 



X?[A] = X?">] [A](X« W ®^- W )[A] -sGf/ [A](2^+ w ®x; w )[A] 



+GSMA] (x« ( %P^- (0) )) [A] + G J [A] (gJL, (xf®v r m ) 8 lf) [A]}. 

(30) 

Two-body terms expansions, obtained from Eq. I|30|l . read as follows: 

X$[A')Xf [A] = X^ <0) [A']Xf <0> [A] 

+ 9 {xf i M<$ [A](xf ® ^ <0) )[A] + Gg,[A'](< ) ^P r <0) )[A']xr (0) [A] 

[A](^+ m ®xf 5 )[A] +Gg,[A'](^ +t °' ^l^IA] 
see appendix 1X1 > 

' (31) 



-99< 

X^[A']Xf[A]=Xf"'[A']A7 A 



s{x/V]^[A](xf ' ®^- (0) )[A] + G<[A](xf 5 ®^- <0) )[A']Xf (0) [A] 
5{^ (0) [A]Gj[A](^+ ( %xf 5 )[A] + G^[A'](^+ (0 ' ®xf , )[A']X? (0> [A] 
see appendix El 
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Obviously, the power expansion of the expectation values can be derived from the quantum average of the preceding 
equations, but not as easily as it seems. Indeed, if one formally writes: 

X2'[A']X?[A] =Y^O{a,b)9 a 9 b 

ah , (32) 

(X? [A']X?[A]) =J2C(*,b)9 a 9 b 

ab 

then C(a,b) is not simply equal to (0(a,b)). Actually, C(a,b) depends on all (0(a',b')) for (a',b') < (a,b), and this 
for two reasons: 



• for a given atom a, the frequency correlation functions {Fp[A']F^[A]) are given by 2nS[A' + A]D pq , where 
D pq depends on the stationary values. But the latter are modified by the second atom and, thus, must also be 
expanded in power of g and g. This implies, for example, that the first term X" m [A']X" ° [A] in the expansion 
of X"[A']X?[A] (Eq. '(EH)) will contribute to all coefficients of (JQ?[A']Xf [A]). 

• the Langevin forces acting on two different atoms are correlated since they both originate from the vacuum 
quantum field. More precisely, their frequency correlation functions depend on their relative distance. This 
dependence is analogous to the correlation function of a speckle pattern (resulting from the random superposition 
of plane waves with the same wavelength but arbitrary directions): 

(if [A']*HA]) = 2nS[A' + A]^r^^>* T^^Xf) 

= ~(g + s)2n6[A' + A]T«';vf q T?r(xPx?) (33 ) 

= -^g + g ynS[A' + A]D^. 

Thus, terms like xf [A'}(xf' } ®2^- <0) )[A] appearing in equation l|31[l will also contribute to higher-order 
coefficients in the power expansion of (X§[A'}X^[A\). One must note that, when R -> 0, Vf q -> §<Vg and 
one recovers the single atom correlation functions given by Eq. (|T7Jl , which emphasizes the consistency of the 
present approach. 

Despite these subtleties, it is nevertheless possible to calculate power expansions of the atomic correlation functions. 
More precisely, in order to emphasize the validity of the present approach, we will compare the results obtain from 
the OB equations and from the Langevin approach. Indeed from the atomic correlation functions, the stationary 
solutions can be calculated by inverse Fourier transform as follows: 

{x?x?,) = T^yll dA ' dA ra A ']*f'[A]}. (34) 

As a specific example, the coefficient proportional to g in the perturbative expansion of (X?, [A']Xf[A]) is given by: 



(jrflAW[A])w = (xrMxfiA])^ + [akt mxf ® pr o, )[A])w 

+ (G%[A'](xf 5 ®^- (0) )[A']Xf C0, [A])(°) 



= G^,[A']G%[A](F^,[A']Ff[A])^ + G-f [A] (xf ) (xf [A']V^ m [A]) 
+ G%[A'](xf ) )(V r m [A']xr[A]) 



(35) 



(0) 



where we have used the fact that terms like (X ai ' X^' )(°) (i.e. zeroth order) actually factorize into (X a )(X^) since 
their fluctuating parts necessarily give rise to higher orders in g and g, see Eq. (1331) . The underlined terms correspond 
to the non-vanishing correlations of the quantum vacuum fluctuations evaluated at the two atom positions. 
Finally, separating elastic and inelastic part, one gets: 

(At[A']A7[A])(s) = (2^) 2 5[A'MA](G^[0](^ 

+ 2n6[A' + A] f-l G^[A-]Gg.[A]^ (0) + G< [A\G% [A']Gj_ fc , [A}D^ {xf) (36) 



+ G^,[A']Gl- k [A']G%[A]Dtf ' (X?, )). 
The corresponding stationary solution then reads: 



<^*f> (9) = g"j mxf'xxrKvtn + G%m(xf)^- W ){ xr) 

+ I J ^A (-l G ^A]G^f + G< [A ]G § 3 , [-A]G^ k , [A]< 0) {X f > (3?) 



+ G%, [-A]G^ k [-A]G%[A]Dtf (xf) 



All quantities above only depend on the stationary values without coupling between the atoms and thus can be 
calculated from the single atom solutions. Furthermore, the integration over A can be performed either numerically or 
analytically by the theorem of residues once the poles of G {i.e. the complex eigenvalues of M) are known. Because of 
causality, they all lie in the lower-half of the complex plane. In practice, we have checked that we effectively recover, 
from the preceding expressions, the results obtained from the full OB equations. In particular, the contribution of 
the correlations of the quantum vacuum fluctuations evaluated at the two atom positions (the underlined term) is 
essential to get the correct results. 

The same kind of expressions can be derived for gg terms, but they are slightly more complicated, since they 
explicitly involve three-body correlation functions, more precisely terms like: 



(o) ) (38) 

[A]' 



which require the calculation of three-points Langevin force correlation functions like: 

(39) 



G< [A]G^.,[A']-L J I dAidA 2 (5[Ai + A 2 - A]G%[A 1 ]G^ k ,[A 2 ] (^Fj,[A']F^[A 1 }F^[A 2 ] 



G<[A}G§ kl [A']l- J! dA 1 dA 2 5[A 1 + A 2 -A]G^[A 1 }G^ k [A 1 }G^JA 2 }^^ 

These correlations functions are non-zero even if they involve an odd number of Langevin forces, emphasizing that 
the statistical properties of the vacuum field fluctuations are far from Gaussian. Nevertheless, the explicit expressions 
of the above quantities can be derived (see appendix [BJ|. They lead to rather complicated and tedious formulae for 
the atomic correlation functions at order gg. From that, we get the corresponding stationary expectations values. 
Again, we have checked that we indeed recover the OB results. 



C. Incorporation of an effective medium 



Finally, and in sharp contrast to optical Bloch equations, it is very easy to adapt all the preceding results to the 
case of propagation in a medium with a frequency-dependent complex susceptibility. Indeed, propagation is controlled 
by the complex amplitude g so that the field radiated by an atom at a distance R and at frequency A will be given 
by: 

n+ [A] = ig Vf q ,V-, [A] exp , (40) 

where ^ + [A] is the (complex) scattering mean- free path satisfying the dilute regime condition /c|^ + [A]| 3> 1. The real 
part of l/^ + [A] describes the exponential attenuation of the field during its propagation in the medium while the 
imaginary part describes the additional dephasing induced by the medium. More complicated formulas, accounting 
for possible variations of £ with position, birefringence effects, or even nonlinearities in propagation, can be derived 
in the same spirit. In all preceding equations, leading to the calculation of the correlation functions, any occurrence 
of the dipole operators must then simply be replaced by: 

VT - VT ^i-Jm) (41) 

while keeping the same " medium- free" coupling constant g. In this way, the present approach can be easily extended 
to the situation where the two atoms are embedded in a medium. In the case of a nonlinear medium, this could lead 
to a self-consistent set of nonlinear equations. 

It is important to stress that accounting for the effective medium is rather straightforward in this frequency-domain 
approach but is a much more difficult task in the temporal-domain approach. Indeed, one basic hypothesis for 
deducing OB equations from the Langevin approach - see section IIII Al - is that the light propagation time between 
the two atoms is much shorter than any typical atomic timcscalc. When this condition is fulfilled, it is possible to 
evaluate expectation values at equal times for both atoms, producing the set of closed OB equations. In the presence 
of a surrounding medium, propagation between the two atoms is affected and this basic assumption may fail. If the 
refraction index of the dilute medium is smoothly varying with frequency, then the corresponding propagation term is 



also smoothly varying with frequency and can be factored out. Thus, except for the exponential attenuation, one may 
recover the OB equations where equal times must be used for atoms 1 and 2. On the contrary, if the propagation term 
has a complicated frequency dependence, the problem cannot be simply reduced to OB equations. It will rather involve 
operators evaluated at the other atom, but at different times, thus leading to a much more complicated structure. 
This difficulty may even take place in a dilute medium with refraction index close to unity. Indeed, the important 
parameter is the time delay induced by the medium, itself related to the derivative of the index of refraction with 
respect to frequency. If the medium is composed of atoms having sharp resonances, the effective group velocity can be 
reduced by several orders of magnitude, consequently increasing by the same amount the propagation time between 
the two atoms. Around the atomic resonance line, the typical propagation time delay induced by the medium over one 
mean free path depends on the laser detuning but is of the order of the atomic timescale for the internal dynamics, 
namely T _1 |3J|. In this case, only the full Langevin treatment developed in this paper can properly account for the 
effect of the average atomic medium. Its practical implementation calls for an investigation on its own and is thus 
postponed to a future paper. We must also note that, if the surrounding medium is composed of the same atoms than 
the scatterers, it is not completely clear that propagation in the medium can be described "classically", i.e. that the 
correlation between the Langevin forces acting on the scatterers and the Langevin forces acting on the medium can 
be safely neglected. 

For the rest of this paper, we will consider two isolated atoms in vacuum. 

IV. MAIN RESULTS 

A. Scattered field correlation functions in the CBS configuration 

In the case of a large number of atoms and for a given configuration, the interference between all possible multiple 
scattering paths gives rise to a speckle pattern. When averaging the intensity scattered off the sample over all possible 
positions of the atoms, one recovers the CBS phenomenon: the intensity radiated in the direction opposite to the 
incident beam is up to twice larger than the background intensity and gradually decreases to the background value 
over an angular range A9 scaling essentially as (fc£) _1 , with I the scattering mean-free path. In the present case, the 
averaging procedure is performed numerically by integrating over the relative positions of the two atoms. As will be 
seen below, the far-field condition kR 3> 1 allows for an a priori selection of the dominant terms contributing to the 
CBS signal. 

The field radiated by the two atoms in the direction n at a distance r ^> R 3> A, in the polarization channel e out 
orthogonal to n (e out • n = 0), is given by: 

fi+rthA] = ~re° ut (pi-[A]e- M i +^1A]e-"- R2 ) (42) 
so that the field correlation function in this channel reads: 

(fi- t [n > A']n+ rt [n,A]> - (^) 2 e° ut e° ut |(^+[A']^-[A]) + (Z>*+[A']X> g 2 -[A]) 

+e lfcn - R (X>2+[A']pi-[A]} +e- ,;fen - R (pi+[A']p2-[A])|. (43) 

The CBS effect occurs when the total phase in the interference terms in the preceding expression becomes indepen- 
dent of the positions of the atom. This phase accumulates during the propagation of the incident laser beam to the 
atoms and during the propagation of the radiated field between the two atoms. The phase factor due to the incoming 
laser beam (a plane wave with wave number = fcnj,) can be explicitly factorized out of the atomic operators as 
follows: 

£,a± =v a± e ±ik L -B. a ( 44 ) 

The other components of X, cf. Eq. (|1()J1 . are populations and not affected by this phase factor. In the single atom 
case, the expectation values of the hereby defined operators are independent of the positions of the atoms. 

Defining <fi = k^ • R and 



9i = ge*+ 92 = ge-*, 



(45) 



the Langevin equations l|29|l then become: 

X Q [A] = X Q<0> [A] + g a G< [A] (x Q ® Z)£-) [A] + g a G a ^ [A] (v^+ ® X Q ) [A], (46) 

In the preceding equation, the Green's functions G are now independent of the position of the atoms, so that the 
phase information due to the incident laser beam is entirely contained in the coefficients g a . 

Frequency correlation functions of the Langevin forces, eq. ■ must also be modified accordingly: 

(ff[A']F?[&]) = -l(gp+g a yir5[A' + A]D^. (47) 

Dropping for simplicity, the tilde notation, the field correlation function 143|) . in the backward direction n = — n^, 
becomes: 

<fi- t [-n L , A']fi + ut [-n L , A]) = (^j {(V 1 ^ [A'^- [A]) + (V 2 + [A']V 2 q - [A]) 

+e- 2l4 '(V 2 +[A'}V 1 q -[A}) + e 2i4, {Vl + [A']V 2 -{A])\. (48) 



The configuration average is then performed in two steps. Since we are working in the limit kR 3> 1, the first one is to 
keep only terms with a total phase independent of kR. In the power expansion with respect to the four parameters g\ , 
g2, gi and §2, this simply amounts to keep terms with even powers of g a g a '- This obviously cancels any <f> dependence. 
More precisely, the field correlation function in the backward direction, beside the trivial zeroth order (in g) term, is 
given by: 

(0- ut hn L ,A']fi o + ut [-n L ,A])( 2 ) = (L} ^^{^[A^IA])^ + (V 2 + [A']V 2 q - [A])^ 

+ {V 2 +[A']V l q -[A])^^ + {V 1 +[A']V 2 q -[A]) i - 9 ^A (49) 



kr 



'i[A',A] +C[A',A] 



The preceding field correlation function still depends on the relative orientation of the atoms through the projector 
7> R , so that, in a second step, an additional average over R must be performed. In the preceding equation, the 
first two terms correspond to the usual "ladder" terms L[A', A] (they are actually independent of the direction of 
observation), whereas the two other terms correspond to the usual "maximally crossed" terms C[A', A]: 



L[A' , A] = ^e° ut e° ut |(^+[A']Pj-[A]}( 9151 ' + {V 2 + [A']V 2 q - [A])^)^ 
C[A', A] = je° ut e° ut |(^+[A']pi-[A])( ffl52 ' + (Vl+[A'}V 2 q -[A})^) 



(50) 



B. CBS enhancement factor 



In the case of linear scatterers, the CBS enhancement factor achieves its maximal value 2 (recall that the CBS 
phenomenon is an incoherent sum of two- wave interference patterns all starting with a bright fringe at exact backscat- 
tering) if the single scattering contribution can be removed from the total signal and provided reciprocity holds. This 
is the case for scatterers with spherical symmetry in the so-called polarization preserving channel h\\ h j3(|. 

In this polarization channel, we have calculated the relevant quantities for an evaluation of the CBS enhancement 
factor when no frequency filtering of the outgoing signal is made. We have thus derived the elastic and inelastic ladder 
terms and the elastic and inelastic crossed terms, together with their corresponding frequency spectra, for different 
values of the on-resonance saturation parameter sq = 2\Ql\ 2 /T 2 . This parameter measures the intensity strength of 
the incident laser beam in units of the natural atomic transition line width T, i.e. its compares the on-resonance 
transition rate induced by the laser to the atomic spontaneous emission rate. For a detuned laser beam, the saturation 
parameter is s(S) and is defined as: 



In the following, different values of the laser detuning have also been considered: 



(a) S = 0, s = s = 0.02 (b) S = 0, s = s = 2.00 

(c) S = ST, s = 2.00, s = 0.02 (d) S = 0, s = s = 50.0 

The ladder and crossed terms Q49J) are separated into their elastic and inelastic parts according to: 

L[A', A] = 2tt(5(A + A') {2tt5(A) L el + L incl (A)} 
C[A', A] = 2tt(5(A + A') {2ir5(A) C e i + C inc i(A)} 



(52) 
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Figure 1: Backscattered light spectrum in the helicity-preserving polarization channel h \\ h. The solid lines represent the 
ladder term (average background intensity value) and the dotted lines represent the crossed (interference) term. For both 
terms, the plotted values corresponds to 7i no i(A)/(C tot +L tot ), see Eq.J^3, where C tot +L tot is the total (elastic plus inelastic) 
intensity scattered in the backward direction. The vertical dashed lines indicate the atomic transition frequency. A corresponds 
to the scattered light angular frequency change with respect to the initial laser angular frequency (A = means thus that 
light is radiated at u>l). Graph (a) corresponds to an on-resonance saturation parameter so = 0.02 and a laser detuning 5 = 
; Graph (6) to (so = 2, S = 0) ; Graph (c) to (so = 2, 6 = 5F) and Graph (d) to (so = 50,5 = 0). At low so, the inelastic 
contribution to the total intensity is small and the ladder intensity is almost equal to the crossed one. For a larger saturation 
parameter, firstly the inelastic contribution becomes comparable to the elastic one and secondly, the crossed term becomes 
smaller than the ladder one. For a nonzero detuning, see graph (c), one clearly observes an asymmetry in the inelastic spectrum, 
which reflects the fact that the scattering cross-section of the atomic transition is maximal for resonant light: the symmetric 
inelastic spectrum emitted by a single atom is filtered out when scattered by the other one. At very large saturation (d), the 
structure of the radiated spectrum becomes rather complicated. 



The corresponding inelastic spectra Li ne i(A) and Ci ne i(A) are displayed in figure^ F° r a sufficiently low saturation 
parameter sq, the inelastic contribution to the total intensity is small and the ladder intensity is almost equal to the 
crossed one (see graph QJi). For larger saturation parameters (see graphs ^ and^i), there are two effects : first, the 



inelastic contribution becomes comparable to the elastic one and second, the crossed term is smaller than the ladder 
one. For a nonzero detuning (see graph^:), one clearly observes an asymmetry in the inelastic spectrum, which reflects 
that the scattering cross-section of the atomic transition is maximal for resonant light (indicated by the vertical dashed 
line): the symmetric inelastic spectrum emitted by a single atom is filtered out when scattered by the other one. We 
also observe that the crossed spectrum is much more reduced than the ladder term, highlighting the non-linear effects 
in the quantum correlations between the two atoms. Finally, for much larger saturation parameters (see graph^?), the 
scattered light almost entirely originates from the inelastic spectrum, like for a single atom. However, contrary to the 
single atom case (for which the scattered intensity reaches a constant value) , the total intensity scattered by the two 
atoms decreases when increasing the incoming intensity. Indeed, since the atomic transitions become fully saturated, 
the nonlinear scattering cross-section of each atom is decreasing, resulting in a smaller total intensity scattered by 
the two atoms compared to the one scattered by a single atom. 

The CBS enhancement factor rj is defined as the peak to background ratio. It thus reads: 



tot 



with: 



c 

^= 1+ 7^ ( 53 ) 

L tot = L c \ + Lj°*! = L c \ + [ — — L inc \(A) 

C° — Coi + Cj ° cl — C i + / 7^rCi nc i(A) 

If the CBS phenomenon is reducible to a two-wave interference, as it is the case here, then the enhancement factor 
rj is simply related to the degree of coherence 7 of the scattered light • If the single scattering contribution can be 
removed from the detected signal, and this is the case in the h || h channel, one has simply 77 = 1 + 7 and consequently 
7 = C tot /L tot . The maximal value for ij is 2, meaning that full coherence 7 = 1 is maintained for the scattered field 
since then C tot = L tot . If all interference effects disappear, meaning C tot = 0, ry reaches its minimal value 1 and 
correspon ding ly coherence is fully lost 7 = 0. Furthermore, one can show that in the h \\ h polarization channel, 
L c \ = Cci 32|. Consequently, as soon as C'"*) < m this channel, the coherence of the scattered light field is 

partially destroyed, since then r\ < 2 and 7 < 1. 



Table I: Ladder (average background) and crossed (interference) terms, see Ea. (|52fl . contributing to the light scattered in 
the backward direction in the helicity-preserving polarization channel h || h. The given values are relative to the incoming 
saturation parameter s. At low so, the inelastic contributions are small and almost equal. Thus C tot ~ L tot and the maximum 
enhancement factor 2 of the linear case is thus recovered, meaning that full coherence 7 = 1 is maintained. At larger so, elastic 
and inelastic terms become comparable. For very large so, the contributions from the elastic terms vanish, like in the single 
atom case. The inelastic contributions are also decreasing, reflecting the fact that the probability for the light to be scattered 
by a saturated atom becomes smaller with increasing saturation. Furthermore, the inelastic crossed term is always smaller 
than the inelastic ladder one. This is a signature of a coherence loss 7 < 1 induced by the quantum vacuum fluctuations. 
However, the ratio Cfnci/^inci does not go to zero as so — > 00 but reaches the limit value 0.096 (for 8 — 0). Also, contrary to the 
single atom case, the properties of the scattered light are not solely determined by the saturation parameter s, but additionally 
depend on the detuning 8, as exemplified by cases (a) and (c) , highlighting the role of the inelastic processes. 
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Our results are summarized in tabled At low saturation parameter s , r\ reaches its maximal value 2 and 7=1. 
This is so because the ladder and crossed inelastic components are almost equal as evidenced in^. Increasing s 
reduces further (7'°^ with respect to L\^ v thus decreasing r\ and 7. In the strongly saturated regime, one thus expects 
7 to decrease. However, there is no reason for the ratio C*°*j/L;°*j to tend to zero as sn — > 00. It rather tends to a 
finite value, which depends on the detuning, in agreement with the results published in Furthermore, keeping s 
fixed and decreasing the saturation parameter s, situation (c), r\ increases, as expected, but to a value which strongly 
depends on sq. In other words, contrary to the single atom case, the properties of the scattered light, are not only 



determined by the saturation parameter s |l9j . Indeed, in both situations (a) and (c), s has the same (small) value, 
but the enhancement factor strongly differs, mainly because the inelastic ladder term has increased. This highlights 
the crucial role of the inelastic processes and of the rather complicated quantum correlations between the two atoms. 

This is not however the full story. Depending on the s and S parameters, a rich variety of situations can be observed, 
with various physical interpretations. These are beyond the scope of this paper, which instead concentrates on the 
basic ingredients of the quantum Langevin approach and will be published elsewhere. 



C. Linear response model 



Some insight on the relative behavior of Ci nc i(A) and Li ne i(A) can be found by comparing the respective formulae 
from which these quantities are extracted: 
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There are twice as many terms contributing to the ladder terms as to the crossed terms. A rather simple explanation 
of this fact is borrowed from the usual linear response theory. Indeed, each atom is exposed to two fields : the 
incoming monochromatic field (angular frequency lol, wave vector k^) and the field scattered by the other atom 
(angular frequency ljl + A, wave vector k p ). In the far-field regime R ;g> A, the incoming field is more intense than 
the scattered field. It thus plays the role of a pump beam with angular Rabi frequency Q l > while the second weaker 
field plays the role of a probe beam with angular Rabi frequency fl p . In this case, the response of each atom is simply 



described by its nonlinear susceptibility ■ More precisely, forgetting about polarization effects, we have: 

(57) 



6V+[A] = e-H^L-kr)-** X++[A ] q+ + e -ik P -R„ x+ _ [A j Q - 



SV-[A] = e ik "- R " x_ + [A] n+ + e '(2ki.-k p ).R Q % __[ A ] Q- . 

where the phases due to the light fields have been explicitly factorized. 

As obviously seen, the two terms x+- an d X-+ generate the forward propagation of the probe whereas the two other 
terms x++ an d X-- can generate an additional field in the direction 2k^ — k p provided phase-matching conditions 
are fulfilled. This corresponds to the usual forward four- wave mixing mechanism (FFWM) 15, 22] . If we now replace 
the probe field by the field radiated by the other atom 0, we get: 



JA] = -L ( e - l (^+2k £ -R„-k £ .R^ x++[A]:D - +e i( fc i J -k i .R,) x+ _ [A] p+j 

"j" L J (58) 

SV^JA] = — { e - l(feit ~ kL - R ' 3) X- + [A] V- +e iV*L-R a +kR-k L -R fl ) x __ [A] 23+}. 

Hence the ladder and crossed contributions are given by (dropping for sake of clarity any frequency dependence): 

« e i ^<^-^- 2kR h ++ X- + 'D-V^ + e^<**-^ X++ X-V-V+ 
+ X +-X- + V+ a V- + e^<*«-WVx + -X-KV + 

(59) 

+ x + -x- + ^+ e i(2ki - (RQ - R3)+2fei?) x + -x-^ + ^ + . 

Averaging these expressions over the positions R a and of the atoms while keeping R 3> A fixed, only terms with 
position-independent phases survive, giving rise to: 

« x+ -x- + Kv- 
lW^x ++ x-~v-v+ + x+ _ X - + v+v-. 

This simple model allows to understand clearly why there are twice more terms in the ladder expression than in 
the crossed one. Fields generated in the FFWM process always interfere constructively in the case of the ladder, 
since they originate from the same atom. Of course, in the preceding explanation, we have discarded polarization 
effects and inelastic processes in the nonlinear susceptibilities. Nevertheless, even if in that case the situation becomes 
more involved, the differences between the ladder and crossed expressions still arise from this local four wave-mixing 

process. For example, in the last line of Eqs. I|55|l and we see that the operator (Gj [A]Xj* <0> eg) ) plays the 

role of a generalized nonlinear susceptibility (actually, the standard ones are recovered from the elastic part of X" ( ). 
Thus we recover the same structure as previously depicted, which leads to similar conclusions. 

Finally, as mentioned above, for large saturation parameters sq, even if in that case the total scattered intensities 
(ladder and crossed) are dominated by the inelastic spectrum, we numerically observe that the enhancement factor 
does not vanish but rather goes to a finite limit 1.096 (for 5 = 0). Field coherence is thus not fully erased, which, at first 
glance, could be surprising since the inelastic spectrum is a noise spectrum at the heart of the temporal decoherence 
of the radiated field. But this only means that both crossed and ladder become vanishingly small relatively to the 
incident intensity. Nevertheless, even if it would be hard to derive it analytically from Eqs. (|55|l and (|56|l . they actually 
decrease at the same rate, resulting in a finite (but small) enhancement factor. 



V. CONCLUSION 



In the case of two atoms, even if the quantum Langevin approach leads to calculations more tedious and involved 
than the direct optical Bloch method, it nevertheless gives rise to an understanding closer to the usual scattering 
approach developed in the linear regime. In this way, one also gets direct information about the inelastic spectrum of 



the radiated light. In particular, it clearly outlines the crucial roles played by the inelastic nonlinear susceptibilities and 
by the quantum correlations of the vacuum fluctuations. Furthermore, since the framework of the quantum Langevin 
approach is set in the frequency domain, frequency-dependent propagation (i.e. frequency-dependent mean- free paths) 
between the atoms can be naturally included. 

The next step would be to adapt the present approach to "macroscopic" configurations (i.e. at least many atoms), 
allowing for a more direct comparison with existing experiments 0- This would provide a better understanding of 
light transport properties in nonlinear atomic media where vacuum fluctuations play a role. In particular, for given 
values of the incident laser intensity and detuning, the nonlinear mean-free path becomes negative in well-defined 
frequency windows. This means that light amplification can be achieved in these frequency windows |35l l38||. The 
atomic media would then constitute a very simple realization of a coherent random laser. 
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Appendix A 



The gg terms in Eq. (|31J) read 
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Appendix B: THREE-BODY CORRELATION FUNCTIONS 



1. Single atom case 

The three-body correlation function for the Langevin force reads: 

C abc [A\ A] = dA^A^Ar + A 2 - A]/^]^] (F£ [A']F£ [Ai]F c a [A 2 ] } , (Bl) 

where /[A] and g[A] are regular functions such that the preceding integral is well defined. Going back to the time 
domain, C a h c [A', A] reads as follows: 

C abc [A', A] = ±-JJ dtdt'e iAt e iA,t ' JJJJ dhdi^dt^h + t 2 - t)S(t 3 + t 4 - t)/(ti) S (t 3 ) (FZ(t')F?(t 2 )F?(t 4 )) . 

(B2) 

Then, from the time correlation properties of the vacuum field, one can show that: 

(F?(t')F b a (t 2 )F?(t 4 )) = 4Tl+T«-5(i' - t 2 ) (x:,(t')XZ(t')F?(U)) 

+ 4T q a +T«-S(t' U) {X:, (t')F h a (t 2 )X?, (U)) (B3) 
+ 4T&T£8(t2 U) (FZ{t')XZ{t 2 )X«,{t 2 )) . 

where the T q± are 15 x 15 matrices defined by [X u T>f\ = ±2T? t X j . 

When taken at the same time, the atomic operators (including the identity 1) define a group entirely characterized 
by the group structure constants e ij k , i.e.: 

k 

k 

so that the preceding equation becomes: 



X i (t)X j (t) = ^i j k Xk(t), (B4) 



(F?(t')F b a (t 2 )F?(t 4 )) = 4T^,T^S(t' t 2 )e a , b , u (X^(t')F c a (t 4 )) 

+ 4T£.T£6tf U) (X:, (t')F h a (t 2 )X?, (t 4 )) (B5) 
+ ^ h +T?-6(t 2 - t 4 )e a , b ? (F«{t')X:{t 2 )) . 

Injecting the preceding relations in C(a, b, c) and going back to the frequency domain, we get: 
C abc [A',A] = 4Tl+T b ye Q , 6 ,"i- J f dA.dA^A, + A 2 - A)/[Ai]#[A 2 ] (X^A' + Ai]F"[A 2 ]} 
+ / dA 33 [A 3 ]/[A-A3]^: M [A' + A3,A-A 3 ] 

+ ±n+T^e a , b ? (f:[A']x:[A\) l-JJ dA l dA 2 5{A 1 + A 2 - A)/[Ai]#[A 2 ] 

= ^ £s ,/i J f dA 1 dA 2 5{A 1 + A 2 - A)/[Ai]g[A 2 ]G" t) [A' + A x ] (F?[A' + Ar]^[A 2 ][) 
+ 4T1+T C " C 7^ J dA 3 .g[A 3 ]/[A-A 3 ]^«[A' + A 3 ,A-A 3 ] (B6) 
+ *r&T£e alb ?G° v [A] (F«[A']F«[A]) l-JJ dA 1 dA 2 S(A 1 + A 2 - A)/[Ai]#[A 2 ] 

= 2nS[A + A>]4T%T b \7e a , b , u D™±- J f dA l dA 2 5{A 1 + A 2 - A)f[A 1 ]g[A 2 ]G a uv [-A 2 ] 
+ 4T1+T C " C 7^ J dA 3 g[A 3 }f[A-A 3 ]D b a f c r[A' + A 3 ,A-A 3 ] 
+ 2n6[A + A']4T«tT c 9 c 7e a , b ,"D-Gl[A]i- J f dA 1 dA 2 5(A l + A 2 - A)/^]^], 



where we have introduced the matrix D^ aa [A', A] defined by: 

D M- [A > )A] = i. JJ dAl dA 2 5(A 1 + A 2 -A')(X?[A 1 ]Fg'[A\XZ[A 2 ]). (B7) 

This matrix is calculated using the same strategy (i.e. going back and forth to the time domain) and one finally 
gets: 

D b ik aaa [A',A] = 2n8[A + A'} {g?J0]L«G£JA']D£* + Gf a [A']G% C [0] D" b 

+*I$T£e b , & v D%± j j dA 1 dA 2 8{A 1 + A 2 - A')GU^i\G a kc [^]G a vu [-^i\ 
+4T^T b - b ,e a , b , v D™^ J f dA.dA^A, + A 2 A')G? a [A 1 ]G^ C [A 2 ]G1[-A 2 ]| 
+ (i- |y dA 3 dA 4 6(A 3 + A 4 - A')G«[A 3 ]G£ C [A 4 ]) 

x // dA^A^Ax + A 2 - A') (X", [Ai]F b Q [A]X"[A 2 ])^ . (B8) 

It may seem that we have taken a loop path and that we are back to square one... However, in the last line of the 
preceding formula, we immediately recognize the matrix D^, Qa [A', A]. Thus, the preceding equation is nothing else 
but a linear system for this matrix. More precisely, f^™ a [A', A] is the solution of the following linear system: 

^fH^A] -Ig^A^D^r [A', A] = 4T a [A',A], (B9) 

with 

/f fe >[A'] = 4T*+T£ -L | y dA 3 dA 4 ,5(A 3 + A 4 - A')G« [A 3 ]G£ C [A 4 ] 
J *,aaa [A , ) A] = 27r<y[A + A 'j {g° [o]L«G^ c [A']D^ + G« [A']G£ C [0]L^J 

+±n+T^e b ,/bZ^ J f dA 1 dA 2 5(A 1 + A 2 - A')G« [A^JA^GU-Ai] 



+4r a + a ,T fe -,6 a , b ,^r^ // rfAi^^Ai + A 2 - A')G« [A 1 ]Gg c [A 2 ]G« u [-A 2 ] 



(BIO) 

In the preceding equations, the Green's function G[A] and the diffusion matrix D aa only depend on the Rabi field 
VLl evaluated at the position of atom a. Thus, for any value of A, numerical values of / and J can be computed, 
allowing for a direct calculation of D^ aa [—A, A}. Furthermore, it is not surprising that the matrix / shows up in 
the linear system. Indeed, the Green's function G[A] governs the time evolution of X through a Fourier transform. 
Thus the time evolution of products of operators X;(i)Xj(t) will be simply governed by the Fourier transform of the 
product of two Green's functions G(t)G(t), which is precisely the convolution product found in /. Finally, from the 
knowledge of the matrix D, we can calculate the value of G a f, c [A', A]: 

G abc [A', A] = 2nS[A + A'} {4T%T£e a , v »D%>±- J j dA 1 dA 2 6(A 1 + A 2 A)/[A 1 ] fl [A 2 ]Gl[-A 2 ] 

+ 4T1+TJ C 7^ JJ dA 1 rfA 2( 5(A 1 + A 2 -A)/[A 1 ] 3 [A 2 ]^ Q [-A 1 ,A 1 ] (Bll) 

+ 4T«tTj c 7e a , b ," J D-Gl[A]i- JJ dA 1 dA 2 5(A 1 + A 2 - A)/[Ai]g[A 2 ]| . 

Of course, we recover the global factor 2n5[A + A'], showing that the time correlation function only depends on the 
time difference t' — t (stationary condition). 

2. Two-atom case 



The calculation of quantities like: 

J_/7' . . , . - . . / „. - .-. 

2tt 



C a f c [A', A] = - / / dA.dA^A, + A 2 - A]/[A!] ff [A 2 ] (^[A']^ [A^ [A 2 ] , (B12) 



follows, more or less, the way described in the preceding section. In particular, it also involves the calculation of a 
matrix T)\^ a [A', A] defined as follows: 

£> ^« ( « [A ' >A] = i. JJdA 1 dA 2 S(A i + A 2 -A')(xt[A 1 ]F b [A}X^[A 2 ]) {9) . (B13) 

The latter is also found to be the solution of a linear system, resembling the preceding one (see Eq. JB9|): 

D^ a(S \A\ A] - ir k a a ,A^}D b a ff aiB) [A\A} = jj^V, A], (B14) 

with 

+ 4T h "+(^ (0) )i- | y dAxrfAa^Ai + A 2 - A')G? [^G^ [A 2 ]G"„[-Ai]i)"" <0> 
+ 4T b V(^ (0) )i- J f dA 1 dA 2 6(A 1 + A 2 - A>)G-« [A^JA^J-A^r^ 
- 2Gj i+ JA'pff i- jj dA 1 dA 2 6(A 1 + A 2 - A')G** [Ax]<5? ro [-A 2 ]lf > [A 2 ])^ 

_2G^ + JA'pf <0) i- II dA 1 dA 2 S(A 1 + A 2 - A')Gl^A 2 }(xf ) [A 1 }xf\-A 1 })^y (B15) 
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